
# clear environment
rm(list=ls())

# make sure working directory is set to folder within which this script is saved
getwd()

#### Lucid 2018 Loss Task Data Summaries ####

# clear environment
rm(list=ls())

# load data
la_data <- read.csv('loss_task_data.csv')

# unusable cases = 573
length(which(is.na(la_data$kwloss)))

# summarize loss aversion
summary(la_data$kwloss)

rm(la_data)

#### Lucid 2018 Lex Task Data Summaries ####

# clear environment
rm(list=ls())

# load data
lex_data <- read.csv('lex_task_data.csv')
attach(lex_data)

## summary stats for all trials
length(lex_rt_check) # total words viewed
summary(lex_rt_check) # summarize reaction times
sd(lex_rt_check, na.rm = T) # sd of reaction times
hist(lex_rt_check, main = "Histogram of Raw Lex Task Reaction Times")
summary(lex_correct_check) # summarize correct responses

## summary stats for correct trials only
length(lex_rt_check[lex_correct_check == 1])
summary(lex_rt_check[lex_correct_check == 1])
sd(lex_rt_check[lex_correct_check == 1], na.rm = T)
hist(lex_rt_check[lex_correct_check == 1],
  main = "Histogram of Raw Lex Task Reaction Times (Correct Only)")

## legal english words only
length(lex_rt_check[lex_latency_check != -99])
summary(lex_rt_check[lex_latency_check != -99])
sd(lex_rt_check[lex_latency_check != -99], na.rm = T)
hist(lex_rt_check[lex_latency_check != -99],
  main = "Histogram of Raw Lex Task Reaction Times (English Words Only)")

## non-words only
length(lex_rt_check[lex_latency_check == -99])
summary(lex_rt_check[lex_latency_check == -99])
sd(lex_rt_check[lex_latency_check == -99], na.rm = T)
hist(lex_rt_check[lex_latency_check == -99],
  main = "Histogram of Raw Lex Task Reaction Times (Non-Words Only)")

## positive words only
length(lex_rt_check[lex_latency_check == 1])
summary(lex_rt_check[lex_latency_check == 1])
sd(lex_rt_check[lex_latency_check == 1], na.rm = T)
hist(lex_rt_check[lex_latency_check == 1],
  main = "Histogram of Raw Lex Task Reaction Times (Positive Words Only)")

## positive, correct words only
length(lex_rt_check[lex_latency_check == 1 & lex_correct_check == 1])
summary(lex_rt_check[lex_latency_check == 1 & lex_correct_check == 1])
sd(lex_rt_check[lex_latency_check == 1 & lex_correct_check == 1], na.rm = T)
hist(lex_rt_check[lex_latency_check == 1 & lex_correct_check == 1],
  main = "Histogram of Raw Lex Task Reaction Times (Positive Correct Words Only)")

## negative words only
length(lex_rt_check[lex_latency_check == -1])
summary(lex_rt_check[lex_latency_check == -1])
sd(lex_rt_check[lex_latency_check == -1], na.rm = T)
hist(lex_rt_check[lex_latency_check == -1],
  main = "Histogram of Raw Lex Task Reaction Times (Negative Words Only)")

## negative, correct words only
length(lex_rt_check[lex_latency_check == -1 & lex_correct_check == 1])
summary(lex_rt_check[lex_latency_check == -1 & lex_correct_check == 1])
sd(lex_rt_check[lex_latency_check == -1 & lex_correct_check == 1], na.rm = T)
hist(lex_rt_check[lex_latency_check == -1 & lex_correct_check == 1],
  main = "Histogram of Raw Lex Task Reaction Times (Negative Correct Words Only)")

detach(lex_data)
rm(lex_data)

#### Lucid 2018 Flank Task Data Summaries ####

# clear environment
rm(list=ls())

# load data
flank_data <- read.csv('flank_task_data.csv')
attach(flank_data)

## summaries for all trials
length(flank_data$flank_rt_check) # total images viewed
summary(flank_data$flank_rt_check) # reaction times
sd(flank_data$flank_rt_check, na.rm = T) # sd of reaction times
hist(flank_data$flank_rt_check, main = "Histogram of Raw Flank Reaction Times")
summary(flank_data$flank_correct_check) # summarize proportion correct responses

## summaries of correct trials
length(flank_rt_check[flank_correct_check == 1]) 
summary(flank_rt_check[flank_correct_check == 1])
sd(flank_rt_check[flank_correct_check == 1], na.rm = T)
hist(flank_rt_check[flank_correct_check == 1],
  main = "Histogram of Raw Flank Reaction Times (Correct)")

## summaries of positive congruent trials
length(flank_rt_check[flank_image_check == 1]) 
summary(flank_rt_check[flank_image_check == 1])
sd(flank_rt_check[flank_image_check == 1], na.rm = T)
hist(flank_rt_check[flank_image_check == 1],
  main = "Histogram of Raw Flank Reaction Times (Pos Cong)")

## summaries of correct positive congruent trials
length(flank_rt_check[flank_correct_check == 1 & flank_image_check == 1]) 
summary(flank_rt_check[flank_correct_check == 1 & flank_image_check == 1])
sd(flank_rt_check[flank_correct_check == 1 & flank_image_check == 1], na.rm = T)
hist(flank_rt_check[flank_correct_check == 1 & flank_image_check == 1],
  main = "Histogram of Raw Flank Reaction Times (Correct Pos Cong)")

## summaries of positive incongruent trials
length(flank_rt_check[flank_image_check == -1]) 
summary(flank_rt_check[flank_image_check == -1])
sd(flank_rt_check[flank_image_check == -1], na.rm = T)
hist(flank_rt_check[flank_image_check == -1],
  main = "Histogram of Raw Flank Reaction Times (Pos Incong)")

## summaries of correct positive incongruent trials
length(flank_rt_check[flank_correct_check == 1 & flank_image_check == -1])
summary(flank_rt_check[flank_correct_check == 1 & flank_image_check == -1])
sd(flank_rt_check[flank_correct_check == 1 & flank_image_check == -1], na.rm = T)
hist(flank_rt_check[flank_correct_check == 1 & flank_image_check == -1],
  main = "Histogram of Raw Flank Reaction Times (Correct Pos Incong)")

## summaries of negative congruent trials
length(flank_rt_check[flank_image_check == 2]) 
summary(flank_rt_check[flank_image_check == 2])
sd(flank_rt_check[flank_image_check == 2], na.rm = T)
hist(flank_rt_check[flank_image_check == 2],
  main = "Histogram of Raw Flank Reaction Times (Neg Cong)")

## summaries of correct negative congruent trials
length(flank_rt_check[flank_correct_check == 1 & flank_image_check == 2]) 
summary(flank_rt_check[flank_correct_check == 1 & flank_image_check == 2])
sd(flank_rt_check[flank_correct_check == 1 & flank_image_check == 2], na.rm = T)
hist(flank_rt_check[flank_correct_check == 1 & flank_image_check == 2],
  main = "Histogram of Raw Flank Reaction Times (Correct Neg Cong)")

## summaries of negative incongruent trials
length(flank_rt_check[flank_image_check == -2])
summary(flank_rt_check[flank_image_check == -2])
sd(flank_rt_check[flank_image_check == -2], na.rm = T)
hist(flank_rt_check[flank_image_check == -2],
  main = "Histogram of Raw Flank Reaction Times (Neg Incong)")

## summaries of correct negative incongruent trials
length(flank_rt_check[flank_correct_check == 1 & flank_image_check == -2]) 
summary(flank_rt_check[flank_correct_check == 1 & flank_image_check == -2])
sd(flank_rt_check[flank_correct_check == 1 & flank_image_check == -2], na.rm = T)
hist(flank_rt_check[flank_correct_check == 1 & flank_image_check == -2],
  main = "Histogram of Raw Flank Reaction Times (Correct Neg Incong)")

detach(flank_data)
rm(flank_data)

#### Flanker Effects Calculations ####

# clear environment
rm(list=ls())

# load data
flank_data <- read.csv('flank_task_data.csv')
attach(flank_data)

## calculate overall flanker effect in milliseconds 
#(1 = pos-cong, -1 = pos-incong, 2 = neg-con, -2 = neg-incong)
# effect = 8.47
(mean(flank_rt_check[flank_image_check == 1 | flank_image_check == 2]) - 
    mean(flank_rt_check[flank_image_check == -1 | flank_image_check == -2])) * 1000

## separate flanker effects for positive and negative targets
# pos = -18.67
(mean(flank_rt_check[flank_image_check == 1]) - mean(flank_rt_check[flank_image_check == -1])) * 1000
# neg = 35.60
(mean(flank_rt_check[flank_image_check == 2]) - mean(flank_rt_check[flank_image_check == -2])) * 1000

# bootstrap for confidence bounds on each of these quantities [note: this may take some time and use substantial memory]
set.seed(625)
boots <- list()
for(i in 1:10000) {
  boots[[i]] <- flank_data[sample(nrow(flank_data), nrow(flank_data), replace = T), ]
}

boot_flank <- c()
boot_pos_flank <- c()
boot_neg_flank <- c()
for(i in 1:length(boots)) {
  boot_flank[i] <- (mean(boots[[i]]$flank_rt_check[boots[[i]]$flank_image_check == 1 | boots[[i]]$flank_image_check == 2]) - 
                      mean(boots[[i]]$flank_rt_check[boots[[i]]$flank_image_check == -1 | boots[[i]]$flank_image_check == -2])) * 1000
  
  boot_pos_flank[i] <- (mean(boots[[i]]$flank_rt_check[boots[[i]]$flank_image_check == 1]) - 
                          mean(boots[[i]]$flank_rt_check[boots[[i]]$flank_image_check == -1])) * 1000
  
  boot_neg_flank[i] <- (mean(boots[[i]]$flank_rt_check[boots[[i]]$flank_image_check == 2]) - 
                          mean(boots[[i]]$flank_rt_check[boots[[i]]$flank_image_check == -2])) * 1000
}

# overall flanker effect
mean(boot_flank) # 8.51
quantile(boot_flank, probs = c(0.025, 0.975)) # -5.41 to 22.54

# positive targets
mean(boot_pos_flank) # -18.55
mean(boot_pos_flank) / sd(boot_pos_flank) # -1.84
quantile(boot_pos_flank, probs = c(0.025, 0.975)) # -38.02 to 1.25

# negative targets
mean(boot_neg_flank) # 35.57
quantile(boot_neg_flank, probs = c(0.025, 0.975)) # 15.51 to 55.33

detach(flank_data)


#### Lexical Decision Task to Construct Figure 4 ####

# clear environment
rm(list=ls())

# load data (for word ordering see Qualtrics logic or word_list_final.xlxs)
lex_data <- read.csv('Negativity+Bias_Summer+2018_October+27,+2021_13.38.csv', stringsAsFactors = F) # start with raw data file
lex_data <- subset(lex_data, lex_data$trump == 3) # remove people who failed attn check
lex_data <- subset(lex_data, as.character(lex_data$task) == "lex") # keep only lex task

# initialize matrices to hold word order and reaction times
word_mat <- matrix(nrow = nrow(lex_data), ncol = 70)
rt_mat <- matrix(nrow = nrow(lex_data), ncol = 70)

# fill in matrices
for(i in 1:nrow(lex_data)) {
  rt_char <- lex_data$rt[i]
  rt_vec <- as.numeric(unlist(strsplit(rt_char, split = ",")))
  rt_vec <- rt_vec[-c(1,2)] # remove two -1's from start
  rt_vec <- append(rt_vec, lex_data$temp_rt[i]) # add final trial
  rt_vec <- rt_vec / 1000 # convert to seconds
  
  rt_mat[i, ] <- c(rt_vec, rep(NA, 70-length(rt_vec)))
  
  word_order_char <- lex_data$word_order_num[i]
  word_order_vec <- as.numeric(unlist(strsplit(word_order_char, split = ",")))
  word_order_vec <- word_order_vec[-1] # remove only one -1 from start

  word_mat[i, ] <- c(word_order_vec, rep(NA, 70-length(word_order_vec)))  
}

# two observations with repeated trials
which(!is.na(word_mat[ ,61])) # obs 222 & 698 got one extra trial
which(duplicated(word_mat[698,]))
which(duplicated(word_mat[222,]))

rt_mat <- rt_mat[ ,1:60]
word_mat <- word_mat[ ,1:60]

rt_means <- c()
for(i in 1:30) {
  rt_nums <- rt_mat[word_mat == i]
  #rt_nums <- ifelse(rt_nums > 5, 5, rt_nums)
  rt_nums <- ifelse(rt_nums > 5 | rt_nums < 0.2, NA, rt_nums)
  rt_means <- append(rt_means, mean(rt_nums, na.rm = T))
}

# word frequencies from the SublexUS database (log10wf measure; see 'word_list_file.xls' for order of words in each category)
word_freq <- c(2.9101,3.1258,1.6435,3.7572,4.1245,2.8938,2.7193,2.9845,1.8261,3.2425, # negative words
               2.7818,2.4609,3.1635,2.6675,3.9569,3.3075,1.8451,4.0796,4.2303,3.2238, # positive words
               1.5563,1.7559,2.1461,1.5315,1.3979,2.749,2.5378,2.0334,2.1004,1.4914) # neutral words

cor(rt_means, word_freq)
word_freq_sq <- word_freq^2
fit <- lm(rt_means ~ word_freq + word_freq_sq)
summary(fit)

setwd("..")
pdf(file="Figures/Figure 4.pdf", height = 2.5, width = 3.46, pointsize=6, family = "Helvetica")
par(mar=c(5,4,1,1))
plot(word_freq, rt_means, pch = 16, xlab = "Word frequency (Lg10WF in SUBTLEX - US)",
  ylab = "Reaction time means (in seconds)", ylim = c(0.00, 1.2), xlim = c(1,5))
x <- sort(word_freq)
y <- fit$fitted.values[order(word_freq)]
lines(x, y) 
text(3.7, 0.40, expression("reaction time ~ word freq + "*"word freq"^"2"))
text(3.7, 0.30, expression("adjusted R"^"2"*" = 0.77"))
dev.off()



